# load libraries
library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.1-10, (SVN revision 622)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-3
# import TEAK CHM
chm <- raster("../NEONdata/D17-California/TEAK/2013/lidar/TEAK_lidarCHM.tif")
# test plot
plot(chm, main = "Plot TEAK CHM using raster package")
# test stretch image
image(log(chm), main = "Plot TEAK CHM image \nJust render pixels that stretch to fill the space")
# plot hist to identify 0 values
hist(chm)
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...):
## 32% of the raster cells were used. 100000 values used.
# set 0 values to NA
chm[chm == 0] <- NA
# look at hist again
hist(chm, xlab = "Tree height (m)") # lots of tall trees
# read in aspect
aspect <- raster("../NEONdata/D17-California/TEAK/2013/lidar/TEAK_lidarAspect.tif")
# plot aspect
plot(aspect, main = "Plot of TEAK aspect data")
Are trees taller on north or south aspects?
First, we will set north-facing slopes to 1 and south-facing slopes to 2. * North-facing slopes are 0-45 and 315-360 degrees, Class = 1 * South-facing slopes are 135-225 degrees, Class = 2
# create array of values
class.m <- c(0, 45, 1,
45, 135, NA,
135, 225, 2,
225, 315, NA,
315, 360, 1)
# shape values into matrix
rcl.m <- matrix(class.m,
ncol = 3,
byrow = TRUE)
# let's look at our matrix
rcl.m
## [,1] [,2] [,3]
## [1,] 0 45 1
## [2,] 45 135 NA
## [3,] 135 225 2
## [4,] 225 315 NA
## [5,] 315 360 1
# reclassify raster using matrix
asp.ns <- reclassify(aspect,
rcl.m)
# plot reclassified raster
plot(asp.ns, main = "North- and south-facing slopes")
Save reclassified matrix as output.
writeRaster(asp.ns,
file = "../outputs/TEAK/TEAK_nsAspect2.tif",
options = "COMPRESS=LZW",
NAflag = -9999)
# compress to keep file size as small as possible
We want to create a raster mask that will compare two rasters only at the selected points of north- and south-facing slopes. This masks out the data we are not interested in.
Important: When using a mask to compare two different rasters, need to pay attention to: Are the two rasters the same resolution? Are the locations matched up or is there georeference data in the raster?
# create a mask from asp.ns
asp.ns
## class : RasterLayer
## dimensions : 577, 543, 313311 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 325963, 326506, 4102905, 4103482 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer
## values : 0, 2 (min, max)
# import NDVI layer
ndvi <- raster("../NEONdata/D17-California/TEAK/2013/spectrometer/veg_index/TEAK_NDVI.tif")
# plot it
plot(ndvi,
main = "NDVI for TEAK site")
# mask data
nFacing.ndvi <- mask(ndvi, asp.ns)
# plot it
plot(nFacing.ndvi,
main = "NDVI for north- and south-facing slopes \nTEAK field site")